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We investigate the single-file dynamics of a tagged particle in a system consisting of N hardcore 
interacting particles (the particles cannot pass each other) which are diffusing in a one-dimensional 
system where the particles have different diffusion constants. For the two particle case an exact 
result for the conditional probability density function (PDF) is obtained for arbitrary initial particle 
positions and all times. The two-particle PDF is used to obtain the tagged particle PDF. For the 
general iV-particle case (N large) we perform stochastic simulations using our new computationally 
efficient stochastic simulation technique based on the Gillespie algorithm. We find that the mean 
square displacement for a tagged particle scales as the square root of time (as for identical particles) 
for long times, with a prefactor which depends on the diffusion constants for the particles; these 
results are in excellent agreement with very recent analytic predictions in the mathematics literature. 



I. INTRODUCTION 

Crowding effects are ubiquitous in cells 1 - - large macro- 
molecules in cells reduce the diffusion rates of particles, 
influence the rates of biochemical reactions and bias the 
formation of protein aggregates^. Furthermore, devices 
used in nanofluidics are becoming smaller; crowding and 
interactions effects between particles are therefore of in- 
creasing importance also in this field. 

An example of a system where crowding is dominant 
is the diffusion of hardcore interacting particles (the par- 
ticles cannot pass each other) in one dimension, so called 
single-file diffusion. For single-filing systems the particle 
order is conserved over time (t) resulting in interesting 
dynamical behavior for a tagged particle, quite differ- 
ent from that of classical diffusion. Examples found in 
nature are ion or water transport through pores in bio- 
logical membranes 3 -, one-dimensional hopping conductiv- 
ity and channeling in zeolites^. Furthermore, in biology 
there are examples where the fact that particles cannot 
overtake one another are of importance: for instance, 
DNA binding proteins diffusing along a DNA chain&£>&. 
Single-file diffusion has also been observed in a number 
of experiments such as in colloidal systems and ring-like 
constructions & 10 ' 1 1 One of the most apparent character- 
istics of single-file diffusion is that the mean square dis- 
placement (MSD) ({xq — xt.o) 2 ) (the brackets denote an 
average over thermal noise and initial positions of non- 
tagged particles, xt is the tagged particle position and 
xq~x) is the initial position of the tagged particle) of a 
tagged particle is proportional to t 1 / 2 for long times in 
an infinite system with a fixed particle concentration; 



the corresponding probability density function (PDF) of 
the tagged particle position is Gaussian. The first study 
showing the t 1 / 2 behavior of the MSD and the fact that 
the PDF is Gaussian is found in Ref. [l2|. Subsequent 
studies include Refs. mHmmillll- The * 1/2 -law and 
Gaussian behavior for long times has proven to be of gen- 
eral validity for identical strongly overdamped particles 
where mutual passage of the particles is excluded, for ar- 
bitrary short-range interactions between particles^ Re- 
cently, a generalized central limit theorem was proved 
for the tagged particle motion^ It is interesting to note 
that a mean square fluctuation that scales as i 1 / 2 also 
occurs for monomer dynamics in a polymer within the 
Rouse modeL 21 ' 22 We point out that anomalous scaling 
of the MSD with time, i.e., that ((xr — xr,o) 2 ) is not pro- 
portional to t, can occur also due to long waiting times 
between particle jump events (when the waiting time dis- 
tribution has a divergent first moment) ] 23 > 24 However, for 
such processes the PDF is not Gaussian; the anomalous 
behavior in single-file systems is not due to long wait- 
ing time densities but rather due to strong correlations 
between particles. 

Although much work has been dedicated to single- 
file diffusion of identical particles, fewer studies has ad- 
dressed the problem of diffusion of hardcore particles 
with different diffusion constants. This type of system 
could be of interest, for instance, for protein diffusion 
along a DNA chain (there is a plethora of DNA binding 
proteins). The single-file system with different diffusion 
constants is illustrated in Fig. [T] The particles each have 
coordinates x = (xi, X2, xn) and initial coordinates 
= (#1,0) #2,0j •••! #iV,o)- D ue t° the hardcore interac- 
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FIG. 1: Cartoon of the problem considered in this study: N 
particles are diffusing in a one-dimensional system. Particle 
j (j = 1, N) has coordinate Xj, initial coordinate Xj,o and 
diffusion constant Dj. The particles cannot overtake, hence 
at all times we have Xj < Xj+i. In our analytic calculation for 
the two-particle case the system size is assumed to be infinite. 
In the stochastic simulations we assume a system of finite 
length L, with reflecting boundary conditions at x = ±L/2. 



tion the particles cannot pass each other, and therefore 
retain their order at all times, i.e., 



L 

— < x% < x 2 < 



< Xn < 



L 

2' 



(1) 



where L is the length of the system (and we assumed the 
ends of the system, at ±L/2, to be reflecting). Particle 
j has diffusion constant Dj (j = 1,..,N). The spatial 
distribution of the particles as a function of time is con- 
tained in the TV-particle conditional PDF V{x, t\xn)\ the 
equations for this quantity were given in Ref. [25| (with 
obvious modifications to account for the different diffu- 
sion constants). We are particularly interested in the 
dynamics of a tagged particle with coordinate xq- with 
initial position xq-,0: which mathematically is obtained 
by integrating P(x, t\xo) over all coordinates and initial 
positions except for xt and xr,o> 25]2G 

To our knowledge, the only studies investigating the 
type of single-file system described above are Refs. [27], 
[28l and [29|. In Ref. the particles were assumed to be 
initially placed at the same position. Also, the 'annealed' 
case, where the diffusion constants were randomized be- 
tween the particles for each new ensemble, was consid- 
ered. In Ref. [28| the hydrodynamic behavior of a two- 
component (two different kinds of particles) single-file 
system with boundary injection and extraction were con- 
sidered. Very recently in the mathematics literature, the 
asymptotic behavior for long times of a tagged particle 
in a single-file system with different diffusion constants 
was obtained for the 'quenched' case (i.e., the diffusion 
constants are the same for each ensemble) for hopping 
dynamics on a latticed 

In this study we extend the results from previous stud- 
ies by (i) analytically solving the problem of diffusion of 
two hardcore interacting particle with different diffusion 
constants in which the initial positions for the two par- 
ticles are arbitrary, and valid for all times. The study 
of diffusion with arbitrary initial conditions is important 
in the field of single-file diffusion since in the derivation 
of the ^/ 2 -law it is assumed that the particles are ini- 
tially randomly distributed, (ii) We introduce a new 
fast stochastic scheme tailored for interacting particle 
systems with different diffusion constants, (iii) For the 
general iV-particle case we verify the asymptotic results 
obtained in Ref. [2^ for long times and large N using our 
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FIG. 2: (top) Cartoon of the problem considered in Sec. [TT] 
two particles are diffusing in a one-dimensional system. Par- 
ticle j (j = 1, 2) has coordinate Xj, initial coordinate Xj,o and 
diffusion constant Dj; in general Di Z?2. The particles 
cannot overtake, hence at all times we have xi < X2- In our 
analytic calculation the system size is assumed to be infinite. 
In the stochastic simulations we assume a system of finite 
length L, with reflecting boundary conditions at ±L/2. (bot- 
tom) Phase space region 1Z (the darker upper area, x\ < X2) 
for two hardcore interacting particles. For the analytic solu- 
tions 1Z extends to ±00. 



new stochastic algorithm, and illustrate the behavior for 
shorter times. 



II. TWO DIFFERENT HARDCORE 
INTERACTING PARTICLES 

We consider a system with two hardcore interacting 
particles diffusing in an infinite one-dimensional system, 
see Fig. [5] (top). 



A. Equations of motion 

The particles each have coordinates x = (x±,X2) and 
initial coordinates xq = (2:1,0, 2^2,0)- The hardcore inter- 
action prevents the particles from passing each other: 



00 < X\ < X2 < 00. 



(2) 



We denote the phase-space region spanned by coordi- 
nates x satisfying Eqs. by 1Z, see Fig. [5] (bottom). 
The temporal behavior of the spatial distribution of the 
particles is contained in the PDF V(x, t\xo) which is gov- 
erned by 



dV(x,t\x ) 
dt 



d 2 



d 2 



dx\ 



' dxi 



D ^+ D ^)nS,t\x ), (3) 
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for x € 7Z [P(x, t\xo) = outside 1Z) and Z?i (-D2) is the 
diffusion constant for particle 1 (particle 2). The initial 
condition is 



V{x, 0|a?o) = <5(xi - a;i, )(5(x2 - 2:2,0) 



(4) 



where 5(z) is the Dirac delta- function. The fact that the 
particles cannot pass each other is described by 

d d 

^ Dl dx~~ D 2Q^)V{x,t\xo)\ Xl ^ X2 = 0. (5) 

The above relation is a no flux condition for the normal 
component of the flux vector across the line x\ = x%, see 
Fig. [2] (bottom): Eq. ([3]) can be written as a continuity 
equation dV/dt = — V • J, where the flux vector is J = 
— {x\D\&P/dxi + x 2 D 2 dV/dx 2 ), and x\ (£2) is a unit 
vector in the x\ (x 2 ) direction. The outward normal to 
the Xi = x 2 interface is (see Fig. [2]) h = (xx — x 2 )/^/2 
which allows us to write Eq. j5| as n ■ J'\ Xl =x 2 = 0. 
This reflecting condition guarantees that the probability 
in the allowed phase-space region TZ is conserved at all 
times as it should. 



B. Solution for two-particle PDF 

In order to solve the equations specified in the previous 
subsection we make the variable transformation: 



X 

Di 

q = 2:2 — 2:1. 
Eqs. ©, © and © then become 
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dV(X,q,t) 



at 

dV(X,q,t) 



= D 



->-Y 
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dX" 1 
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V{X,q,t) 



dq 



q=0 



P(X,q,t^0) = jS(X - X )6(q - q ) 



(7) 



w here 7 = [D l + D 2 ) / '(2x/Dx~D^), X = y/D 2 /DxXx,o + 
y/ D\/ D 2 X2fi and go = 2:2.0 — 2:1,0 and we introduced the 
effective diffusion constants 



D x = 



D x +D 2 



D q = D x +D 2 



(8) 



For the case of identical diffusion constants, D\ = D 2 = 
D the equations above express the fact that the rela- 
tive coordinate q diffuses with a diffusion constant 2D, 
whereas the center-of-mass coordinate X diffuses with a 
diffusion constant D/2^ 

Eq. ||7J) allows a product solution of the form 



V(X,q,t)=P Ji (X,t)V ! (q,t) 



(9) 



where 



V x {X,t) 



(AirD x ty/ 2 



exp(- 



(X - x f 

4D x t 



(10) 



and the solution for r P q (q,t) is obtained via the method 
of images^ according to 



V(q,t) = 9(q) 
x I exp(— 



1 



(AttDH) 1 / 2 

(q - <7o) 2 x . 



exp(- 



(g + q f 
ADH 



(11) 



where 6(q) is the Heaviside step function, 9{q > 0) = 1 
and 0{q < 0) = 0. Returning to our original coordinates, 
Eq. (J9]), (fT0|) and (fTTjl become, after some algebraic ma- 
nipulations: 



V{x,t\xo) = ®{X2 - xi) 



1 



{AnDxt) 1 / 2 (4 7 rD 2 t) 1 / 2 



x [exp(- 



(x 1 



x i,o) 



■) exp(- 



(x 2 - x 2fi f 



ADnt 



I {Xx-X\ fi f (x 2 -4. ) 2 

ex P( TFTT ) ex P( 77T1 )I 12 ) 



AD x t ' v AD 2 t 
where the effective image initial positions are 
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2D, 



Dx 



D 2 
D 2 



Dx + D 2 



X2M 



X2,0 



(13) 



Notice that we have x\ — 2^ = —(2:2,0 — £1,0)1 i-e. the 
distance between the image initial positions is the same 
as the distance between the initial positions. Eq. (1131) is 
a non-trivial extension of the image positions for identical 
particles or for a system where the particles initially start 
out at the same point in space: For Dx — D 2 we have 



x 2 and x 



2.0 



X1.0 as it should^ For the case 
2:1,0 = x 2,q = the results above reduce to the results 
obtained in Ref. Note that, in contrast, when D\ ^ 
D2 and 2:^0 7^ 2:2,0 the image initial positions, Eq. (|13p . 
depend on Dx and D^. 

It is interesting to compare the above result for 
V(x,t\xo) to that of a Bethe-ansat a 32 ' 33 . It is straight- 
forward to show that Eq. (fT2"| can be written: 



T(x,t\x ) = 9{x 2 



dkx r°° dk 2 



-lk 2 X2 



x |e" ClXl e lK2X2 + g(k 1 ,k 2 ,x 1 , x 2 )e lk2Xl e lklX2 } 



(14) 



where 
g(k 1 ,k 2 ,x 1 ,x 2 ) 



cxp| 



r Di-D2 
Dx+D 2 



(kx + k 2 )(x 2 - Xl )} (15) 



We note that Eq. ([T4]) has the form of a Bethe ansata 3 ^, 
where the "scattering coefficient" g depends on xx and 2:2 
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(in the standard Bcthe ansatz the scattering coefficient 
only depends on ki and A^); the standard Bcthe-ansatz 
satisfies the equations of motion and the boundary con- 
ditions for fixed k\ and in contrast, the solution above 
does not - it is only after the integrations over k\ and k% 
are performed [with the appropriate x\ and xi dependent 
"mixing" of k\ and ki from the 2nd term in Eq. (|14[) ] 
that the correct solution for V(x,t\xo) is obtained. For 
the case of identical diffusion constants D\ = D2 = D 
the mixing of k\ and fc 2 in g is absent and we have g = 1 
in agreement with previous studies^ 



C. Tagged particle PDF 

By integrating the two-particle PDF we obtain the 
tagged particle PDF. The tagged PDF (for fixed 
initial positions) for particle 1 is pi(xi,t\xo) = 
dx2'P(x, t\x*o). Explicitly using Eq. (fT2")h we have: 



pi{xi,t\x ) 



exp(- 



(xi - x lfi f 



(4 7 r J D 1 i) 1 /2 

cxp( ~ 

x-erfc( — 1 ) 
2 y/ZDrf 



W\t 



(xi - x\ fi f 
W x t 



(16) 



where erfc(z) = 1 — erf(z) is the complementary er- 
ror function, with erf(z) = (2/y^vr) J Q G?yexp(— y 2 ) be- 
ing the error functional and x\ and x 2 are given 
in Eq. (fHi|) . The tagged particle PDF for particle 2 
Pi{xi-,t\x<i) = Jl^dxiViXjtlxo) is obtained by the re- 
placements xi *-* —X2, x\fi <-> — X2,q, i ! 10 <-* -^2,01 an d 
£>i «-> D 2 in Eq. (16]). We point out that Eq. (HJ) does 
not give an MSD which scales with time as t 1 / 2 (see In- 
troduction); it is only in the limit of a large number of 
hardcore interacting particles that Pj(xj, t\xo) (with an 
additional average over the initial position of non-tagged 
particles), for a center particle, becomes a Gaussian with 
a width that scales as the square root of time, see next 
section. 

A limit not accessible through previous approache a 27 i 31 
is that where one of the particles is immobile D 2 = 0: set- 
ting x 2 o = for convenience and taking the limit D2 — > 
Eq. (fTo| becomes 



pi(x 1 ,t\x )\ D2 =o = 0(-xi) 



1 



(4ttL>i0 1 /2 



x I ex P("^ TFTL ) + CX P( )l ( 17 ) 



4£>ii 



W^ 



in agreement with the diffusion of a particle near a re- 
flecting wall as it should. We have above used the fact 
that erf(±oo) = ±1. 
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FIG. 3: Tagged particle PDF pj(xj,t\xo) for two hardcore 
interacting particles (j — 1,2). The solid and dashed left- 
most (blue) [rightmost (red)] curves corresponds to the tagged 
particle PDF for particle 1 [particle 2] at different times as 
given in Eq. (|16[) . The ratio of the diffusion constants is 
D2/D1 — 0.3. The symbols correspond to the results of the 
Gillespie simulation, ensemble-averaged over n ons = 50000 en- 
sembles for M — 500 lattice sites; the data was binned into 25 
bins. The vertical bars at the bottom of the figure corresponds 
to the initial positions for the two particles (xi,o = 0.10L and 
£2,0 = 0.15L). For long time t > L 2 /D the equilibrium is 
reached - the dashed lines correspond to the analytic result 
as given in Eq. (|18)>. 



In Fig. [3]wc illustrate the results for the tagged particle 
PDFs as given in Eq. (| 16|) . We compare to stochastic 
simulations using a new stochastic algorithm described 
in Appendix [X] In the simulation we assume a finite 
box; the main effect of the finite box (with reflecting 
conditions) is to modify the long-time limit (i.e. for t 3> 
L 2 /Dj), to the following equilibrium PDFs>^ 



(18) 



The results given in Eq. (fT8"|) is obtained by di- 
rect integration of the two-particle equilibrium PDF 
P eq (xi,a;2) = 29(x2 — xi)/L 2 . The equilibrium results 
are independent on D\ and D 2 as it should. In Fig. [3] 
we illustrate the result of a Gillespie simulation using 
Hens = 50000 ensembles on a lattice with M = 500 lat- 
tice points, and compare to the PDF Eq. fTo]) as well 
as the equilibrium PDF, Eq. (|18[) . We notice excellent 
agreement within the limit of applicability. 

It remains a challenge to generalize the results in this 
section to the ./V-particle case and arbitrary times. In 
the next section we perform stochastic simulations for N 
particles and verify the asymptotic results in Ref. [2^ for 
long times and N large. 
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III. N DIFFERENT HARDCORE 
INTERACTING PARTICLES 

For N identical point particles we have the standard 
result ((x T - x T ,o) 2 ) = (1/ Q)(4Dt/n) 1 / 2 for the MSD of 
a tagged particle, where g = N/L is the concentration 
of particles (N, L — > oo with g kept fixed) , 12 i 13 i 25 For 
single-file particles with different diffusion constants very 
recent results show that the motion for a tagged particle 
stochastically jumping on a lattice (exponential waiting 
time between jumps) is a fractional Brownian motion (so 
that the PDF is Gaussian)^ The MSD was in Ref. [H 
shown to take the same form as for identical diffusion 
constant but where D above is replaced by an effective 
diffusion constant, i.e., we have: 



((x r - xt,o) 2 ) = « 



^D eS t 



1/2 



with the prefactor 



1-/ 



(19) 



(20) 



where / = N/M, M is the number of lattice points and 
a the lattice spacing (N,M — > oo with / fixed). In the 
continuum limit (lattice spacing a — > 0) we get the point- 
particle result n = l/g. In the continuum limit but 
with finite-sized particles we have, as in Ref. I2EL that 
K = (1 — gA)/g, where A is the size of the particles. 
The effective diffusion constant appearing in Eq. (fT9|) is 
obtained by averaging the friction coefficients (inverse of 
diffusion constants) according to: 



1 1 N 1 



(21) 
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FIG. 4: Mean square displacement for a tagged particle in the 
center of a single-file system with different diffusion constants. 
The upper (blue online) marks show the simulation results 
for a system where the diffusion constants are 'quenched' and 
random with 0.01D ma x < Di < £> max (the inset shows the 
diffusion constants used). The lower (red online) marks are 
the results for an alternating sequence of diffusion constants: 
Di = D m ax,I>2 = 0.01D max ,D 3 = -Dmax, D 4 = 0.01D ma x etc. 
For both cases the tagged particle's diffusion constant Dr was 
set equal to D m ax- The solid blue and dashed red line are the 
analytic result as given in Eqs. (|19|l - (|21| l. In all simulations 
the tagged particle was initially placed at the center lattice 
point and the non-tagged particles then randomly positioned 
to the left and right of the tagged particle. The errorbars 
are the standard errors. The following parameters were used: 
M = 10001, N = 2001 and the number of ensembles n cns = 
3200. Without loss of generality, we set a — 1 and Dr = 1 
(i.e., a and Dr determine the units of length and time in the 
problem) in all simulations. Particle 1001 was taken to be 
tagged. 



provided the limit on the right-hand side exists^ 5 - The re- 
sult above is obtained for the (realistic) 'quenched' case, 
i.e., Di are kept the same for all the ensembles. The 
results above are thus stronger than the results in Ref. 
I27I where the 'annealed' case, i.e. for each ensemble the 
diffusion constants are reshuffled between the particles, 
was studied. We also point out that the results above 
are valid for an initial equilibrium density of particles, 
whereas the results in Ref. [2?] are limited to the case 
that all particles are initially placed at the same point in 
space. 

In Fig. 2] we show results of stochastic simulations for 
TV = 2001 particles, with the middle particle (particle 
number 1001) being tagged. The tagged particle is ini- 
tially placed at the center lattice point and the remaining 
particles are randomly positioned (avoiding multiply oc- 
cupied lattice sites^) to the left and right of the tagged 
particle for each ensemble. The details of our stochastic 
scheme is presented in Appendix [X] Two cases are pre- 
sented: the upper (blue) marks shows simulations for the 
case of random 'quenched' distribution of diffusion con- 

Thc tagged 



particle diffusion constant Dj- was set to D maxi and the 
diffusion constants used are shown in the inset in the fig- 
ure. The lower (red) marks represent results for an alter- 
nating set of diffusion constants: the first particle has dif- 
fusion constant -D m ax, the second particles has 0.01Z? max 
the third -D max etc. The tagged particle has diffusion 
constant D max . For short times, t <C l/(g 2 Dr), we see 
that the case of 'quenched' random and alternating dif- 
fusion constants give the same MSD (the tagged particle 
has a diffusion constant equal to -D max in both cases). 
There has been few collisions between particles and the 
MSD is proportional to t, see Fig. HI The MSD in the 
short-time regime is slightly smaller than that of a free 
particle ({xt — xt.o) 2 ) = 2Dq-t [dotted line], simply due 
to the fact that in some ensembles the tagged particle will 
initially have a non-tagged particle at a neighbouring lat- 
tice site (every fifth lattice site will on average contain a 
particle in the simulations in the figure). For long times, 
t ^> l/(g 2 D e g), there is a cross-over to a single-file regime 
with the MSD proportional to t 1 / 2 . We notice an excel- 
lent agreement with the stochastic simulations and the 



prediction in Eq. (TT9l) - ([2T1) [solid blue and dashed red 
line] for long times. We point out that the average diffu- 
sion constants [(1/iV) Yli=i A ] for the two cases above 
are very close (more precisely, the two cases converge to 
the same average diffusion constant for N — > oo); Fig. 
2] thus clearly illustrates that it is the average friction 
coefficient which determine the long-time behavior for 
the system rather than the average diffusion constant. 
For very long times (beyond the time window in Fig. 0]) 
and finite L, the equilibrium PDF for the tagged particle 
should be reached, see Ref. [25| for an explicit expression. 



IV. SUMMARY AND OUTLOOK 
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FIG. 5: Schematic illustration of the lattice on which N par- 
ticles jump in the Gillespie simulation. The number of lattice 
sites are denoted by M and the lattice spacing is a. Parti- 
cle i jumps to the left (right) with rate fei-2 (fei-i). Note 
that these rates change during the simulation: if particle 1 is 
at the leftmost site then we impose the reflecting condition 
fei = 0. Similarly if particle N is at the rightmost site we have 
k2N-i = 0. If particle i and i + 1 are at neighboring sites we 
set fc2i-i = &2i = due to the hardcore repulsion. For the 
remaining configurations we have that the hop rates take the 
"free" values: fc M = kl. 



In this study we have investigated the (single-file) dy- 
namics of hardcore interacting particles with different dif- 
fusion constants diffusing in a one-dimensional system. 
For the two particle case we obtained an analytic re- 
sult for the conditional PDF (for arbitrary initial par- 
ticle positions and all times), from which we calculated 
the tagged particle PDF, see Eq. (|f 6[) . For the general 
./V-particle case an asymptotic expression for the mean 
square displacement of a tagged particle for long times 
was given, Eq. (flU)) , and excellent agreement was found 
with our new computationally efficient stochastic simu- 
lation technique based on the Gillespie algorithm. 

It will be interesting to see whether it is possible to 
generalize our two-particle PDF to N particles with dif- 
ferent diffusion constants, in order to access the full time 
behavior, and thus going beyond the asymptotic results 
in Ref. [2!J We point out that the TV particle results given 
in this study assumed the mean friction constant to be 
finite; we are currently considering the case of a distri- 
bution of friction constants with diverging first moment. 

The problem studied here, the dynamics of interact- 
ing species of different kinds, shares many features with 
the dynamical behaviour of cellular (and other biolog- 
ical) systems where heterogeneity and interactions are 
important factors. We hope that our study will inspire 
to further studies of many-body biology effects in living 
systems. 
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APPENDIX A: NEW EFFICIENT STOCHASTIC 
ALGORITHM BASED ON THE GILLESPIE 
ALGORITHM 

Stochastic simulations using the Gillespie 
algorith m 37 i 38 i 39 is a convenient technique for gen- 
erating stochastic trajectories for interacting particles. 
Briefly, we consider hopping of N particles on a lattice 
with M lattice sites. The dynamics is governed by the 
'reaction' probability density (rPDF) 

2JV-1 

P(t,h) = fc Al exp(- k » T ) ( A1 ) 

where r is the waiting time between jump events and 
fc M are the jump rates. There are 27V jump rates for 
the process considered here (/i = 0, 2N — I): the rate 
for particle 1 to TV jumping to the left and right respec- 
tively. We enumerate these rates such that k 2 i-2 is the 
rate for particle i jumping to the left, and fei-i is the 
rate for particle i jumping to the right, see Fig. [5] For 
the case that a particle have no neighbors nor are at the 
end lattice points we set fc M = where kj[ are the "free" 
hop rates. For the case that particle I (particle N) is 
at the leftmost (rightmost) lattice point we have the re- 
flecting condition k\ = (k 2 N-i = 0). If two particles 
are at neighboring sites the hardcore repulsion requires 
the right (left) rate for the leftmost (rightmost) particle 
equals zero, i.e. if particle i and i + 1 are at neighboring 
site we set k 2 i-i = k 2 i = 0. Jump rates that are not 
set to zero are equal to their "free" hop rates. Thus, a 
stochastic time series is generated through the steps: (I) 
place the particles at their initial positions; (2) From the 
rPDF given in Eq. (|A1[) we draw the random numbers 
r (waiting time) and /x (which particle to move and in 
what direction); (3) Update the position of the chosen 
particle, the time t and the rates kj for the new configu- 
ration and return to (2); (4) The loop (l)-(3) is repeated 
until t > t s top, where t s top is the stop time for the simu- 
lation. This procedure produces a stochastic time series 
for < t < tstop- If steps (l)-(4) are repeated n cns times 
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one obtains a histogram of particle positions (at specified 
times). The ensemble averaged results of a Gillespie time 
series is equivalent to the solution of a master equation 
incorporating the rates given above) 37 i 38 see for instance 
Refs. Ho| or [32| for the explicit expression for this mas- 
ter equation. In the limit a — > with fixed diffusion 
constants Di = (fefi-2 + ^2i-i) a2 /2 the master equation 
approaches the diffusion equation (assuming no drift, i.e. 



that k 



f 

2i-2 



k 



2j-i) as specified in section |TT] for an infi- 



nite system. 

In the two subsequent subsections we consider two dif- 
ferent methods for generating a stochastic time series for 
the type of dynamics described above: (A) the direct 
method and, and our new approach (B) the trial-and- 
error method. We find that the latter method is superior 
in computational speed whenever the fraction of occupied 
lattice sites is not close to one. 



2. (B) The trial-and-error method 

In this subsection we present a new method for gener- 
ating a time series according to the rPDF as given in Eq. 
(IAI I) ; we call this method the trial-and-error method. 

1. Generate the partial sums of the free rate constants 



Po = 

Pi* = J2 k -' M=l,-,2iV 



(A4) 



2. Generate an initial configuration of particle posi- 
tions. 

3. Set the waiting time equal to zero, r = 0. 



1. (A) The direct method 

Let us now review the "direct method" used for ob- 
taining a stochastic time series using the rPDF as given 
in Eq. (|Aip . Briefly, the direct method as applied to the 
present type of system involves the following step a 37 ' 38 i 39 

1. Place the N particles at their initial positions and 
assign free hop rates kl. 

2. Draw two random numbers r\ and ri (0 < n, T2 < 
1) 

3. The waiting time is obtained from the first random 
number as: 



(A2) 



4. The 'reaction' /i is determined by using ri to de- 
termine the [i which satisfies 



E 

i/=0 



2N-1 



K < r 2 ^ K < 



v=0 



(A3) 



5. Use the obtained /x and r to update the particle 
positions Xi(t) and the time. After the chosen par- 
ticle has been moved, check the local environment 
for the particle and if necessary update the rate 
constants k^ accordingly {k^ is either equal to kj L 
or 0). 

6. Return to step (2). 

This procedure produces a stochastic time series Xi(t) 
for the particle positions. If steps (l)-(6) are repeated 
n cns times (n cns ensembles) one obtains a histogram of 
particle positions (at specified times). 

The direct method as applied to the present type of 
system is computationally slow, since at each step the 
sum of all rate constants has to be recalculated. 



4. Draw a random number n (0 < r\ < 1). The 
waiting time is updated according to: 



1 , f 1 
t H log — 

P2N V r l 



(A5) 



5. Draw a random number r 2 (0 < T2 < 1). A trial 
'reaction' is determined by the /i which satisfies 



Pii < r 2 p2N < Pii+i 



(A6) 



6. Check if there is a reflecting wall or neighboring 
particle preventing the 'reaction' /i from occurring, 
if so return to step (4). If fi is allowed then update 
the particle positions and the time. 

7. Return to step (3). 

The scheme above produces a stochastic time series 
Xi(t), and if steps (2)-(7) are repeated n ons times (n ons 
ensembles) one obtains a histogram of particle positions; 
note that step (1) does not have to repeated for each en- 
semble. An efficient method for performing the search for 
the fi satisfying the inequality in Eq. (|A6|) is presented in 
Appendix [Bj Notice that in the scheme above one does 
not have to update the k^ at each time step, since we are 
only using the free rate constants kj^ (which are fixed at 
all times) in the present scheme, nor do we need to per- 
form a sum over all rate constant at each time step. The 
only price we have to pay compared to the direct method 
is the rejection step (6), which may cause us to repeat 
steps (4) and (5) several times. However, whenever the 
fraction of occupied lattice points is not close to one (or 
more precisely, essentially whenever the occupation frac- 
tion is such that we can find an allowed reaction in less 
than N steps) we expect that the trial-and-error method 
is faster than the direct method. A proof that the di- 
rect and the trial-and-error methods are mathematically 
equivalent is presented in Appendix [Cl 
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APPENDIX B: FAST METHOD FOR 
DETERMINING /x FROM EQ. (fA6|l 



(|A1|) . Let us first define a trial-and-error relaxation time 
as 



In this appendix we give a simple, yet fast, algorithm 
for determining the fi satisfying Eq. (|A6[) . We start by 
noticing that if all are identical then we can directly 
satisfy Eq. (|A6j) by choosing fx = [2iVr2], where [z] gives 
the largest integer smaller than z. The algorithm below 
uses this fact to directly find the correct [i for identical 
rate constants; for non-identical rate constants we expect 
the number of attempts before the correct fj, is found to 
scale not worse than logiV (see below). The algorithm 
is: 

1. Initialize the two "boundary" parameters /j} elt and 
^ right to p) clt = and ^ right = 2N. Also, introduce 
p and p rlght which are initialized to the largest 
and smallest of the partial sums, i.e. we set initially 
p left = Pa = and p right = p 2N , see Eq. (|A4)) . 

2. Make a guess for fi using the formula 



tte 



left. 



< r 2 p 2N -P 

L fright picft 



(^"-O + Mhft]. (Bl) 



m, if 



3. Check if the guess value for /x satisfies Eq. 
so terminate the loop. 



4. If the guess value for fi does not satisfy Eq. (1A0I) 
we separate between the cases: (a) if r 2 p 2 N < 
M then move the right boundary parameters 
by changing // lght -» /j guoss and p rlght -> p Mguess ; 
(b) if r 2 p 2 N > TVguoas+i then move the left bound- 
ary parameters by changing /i lcft — > ^mss + 1 and 



P 



left 



/^gucss 

-i. Then return to step (2). 



The algorithm above will narrow down the search to a 
smaller and smaller segment along the /x-axis, until we 
manage to obtain the correct fi. The initial guess will 
always be /x guess = [2Ar 2 ], see Eq. (|B1I) . and hence for 



v 2iv-i uf 



(CI) 



and we note that for each (successful or unsuccessful) 
attempt we draw a waiting time from the PDF: 



Pte(t) = -^—e 
tte 



-t/tte 



(C2) 



see step (4) in the trial-and-error Gillespie scheme. To 
construct P(t, (i) within the trial-and-error method we 
have to wait for a successful attempt. We therefore sep- 
arate between the following (mutually exclusive) events: 

• success on the first attempt: the probability for this 



= Q = , t = — (C3) 

TG 



where tq — l/E^=o 1 ^m] 1S relaxation time for 
a successful event. 

• one failed attempt, then success: the probability 
for this is o\ = q(l — q). 



• m failed attempts and then success: the probability 
is o m = q(l - q) m . 



Given one of the above events the probability that reac- 
tion /.t occurs is 



hi. 



2^=0 K » 



(C4) 



identical free rate constants k f = k[ = ... = A^_ a the The successful rPDF thug becomes 



algorithm above directly finds the correct /i. To argue for 
the maximally log A-scaling for the number of iterations 
note that if we had chosen /x guess = [(// Ight — /i )/2] in 
the algorithm then the number of steps to find the correct 
H would be maximally [1 + log 2 N] . The algorithm above 
should perform better than this because the values of 
/ig UCSS should approach the correct \i faster. 



APPENDIX C: PROOF OF EQUIVALENCE OF 
THE DIRECT METHOD AND THE 
TRIAL-AND-ERROR METHOD 

Let us finally prove that indeed the trial-and-error 
method is equivalent to the Gillespie rPDF, Eq. (|A1[) . 
More precisely we prove that the successful 'reaction' 
probability density function P(r, fi) agrees with Eq. 



P(t,h) = p m ct pte(t) 

+P^Oi I drip T E(Tl)/OTE(T - Tl) 
JO 

+P(iO- 2 J dTldT 2 p TE (Tl) 

xPte(^2 - ti)pte(t - t 2 ) + .. (C5) 

The first term is the PDF for a successful move of type 
fj, on the first attempt. The second term is the PDF for 
a successful move of type /i on the second attempt, and 
is obtained by considering the waiting time density to 
first make an unsuccessful attempt at time t\ followed 
by a successful attempt at time r; since T\ can take any 
value between and r we must integrate over this inter- 
val. Similar arguments give expressions for higher order 
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terms. In order to express the result given in Eq. (|C5P in is a geometric series E^°=o a " = V(l — a ); the explicit 
closed form wc make a Laplace-transform which gives: 



due UT P(r,p) = p^aopTniu) 
+P ai cti[pte(w)] 2 + P m ct 2 [pte(w)] 3 + 

OO 

n=l 

oo 

Pt,qpTE{u) ^[(1 - q)pTE{u)] n 



(C6) 



where pte(w) = 1/(1 + utte) is the Laplace-transform 
of Eq. (|U2j) . Using the fact that the series in Eq. (|C6)) 



expression for pte( u ) and Eq. (|C3p . we finally find: 

fc, 



P(ti,//) = 



1 + UTTE/q 



(C7) 



In the time domain this is the same as Eq. (|Alj) . which 
thereby completes the proof. 
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